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In the paper a bridge functional for the closure relation of the three-dimensional reference interaction model 
(3DRISM) equations is proposed. The effectiveness of the bridge for the aqueous solutions of the carbon 
nanomaterials is tested. In the paper two classes of systems are investigated: (i) infinitely diluted aqueous 
solutions of the Lennard-Jones (LJ) spheres and (ii) infinitely diluted aqueous solution of the carbon nan- 
otubes(CNT). The bridge functional is fitted to the molecular simulation data. It is shown that for all the investi- 
gated systems the bridge functional can be approximated by the exponential function which depends only on 
the solute/solvent size ratio. It is shown that by using the proposed bridge functional is possible (i) to predict 
accurately the position of the first peak of the water oxygen density distribution function (DDF) around the 
solute, (ii) to improve the accuracy of the predictions of the first peak's height of the water oxygen DDF around 
the solute (iii), to predict correctly the water hydrogen DDF behavior in the vicinity of the CNT. 
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1. Introduction 

Integral equation theory of liquids (IETL) allows to calculate the thermodynamic parameters 
of a liquid and an amorphous matter by using the local microscopic structure which is deter- 
mined by the density distribution functions (DDFs) DDFs are straightforwardly connected to 
the total correlation functions which in turn can be calculated from the Ornstein-Zernike (OZ) 
equation ant the closure relation 0, For the molecular systems the correlation functions in 
a general case depend on the six variables. Solving the six-dimensional molecular OZ equation 
is still a challenging task. In practice usually the approximations of the OZ equation are used 
such as the reference interaction site model (RISM) |3| and the three-dimensional reference in- 
teraction site model (3DRISM) 0]. For RISM and 3DRISM models the effective computational 
algorithms were developed [5 -9]. The RISM and 3DRISM equations allow to predict the properties 
of the liquid in a thermodynamic limit. In many applications they can be used as a substitution of 
the time-consuming molecular simulations. The 3DRISM was recently implemented in the Amber 
molecular modeling package and in the Amsterdam Density functional (ADF) quantum mechanics 
softwarejljj, U|. This fact reflects the constantly increasing interest to this method in the chemical 



community. The 3DRISM can be used for investigating the structure of the molecular complexes 
[l2 - 14 1 , for determining solvation effects on conformational equilibria [TBI, [l(| , for the accurate 
predicting of the solvation free energy (l7l - l23l | , and in many other applications [U 27 1 . Although 



the RISM and 3DRISM equations predict qualitatively correctly the solvent structure. They often 
fail to reproduce the qualitative parameters such as the radial distribution functions (RDFs). The 
main source of errors is not the model's approximations but the closure relation. The rigorous clo- 
sure relation includes the so-called bridge functional which have in its analytical representation the 
infinite series of the integrals with the growing dimensionality 0. This makes the bridge functional 
practically incomputable. In practice the approximate closures are used. This causes the inaccu- 
racies in the calculations. There was proposed a number of closure approximations. To name a 
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few most popular: the Hyper Netted Chain (HNC) closure where the bridge functional is omitted, 
the Percus-Yevick (PY) closure for hard spheres [28j, Martynov-Sarkisov closure [29j, the Verlet 
Modifiet closure [3(j], and others[31]. The most of the functionals describe more or less accurately 
the simple liquids of hard spheres or LJ spheres. But even for the two-component solutions of LJ 
spheres of different size the effectiveness of the standard bridge functionals is questionable. 

Another way to obtain the universal bridge functional is a parameterizing of the bridges which 



are obtained from the molecular simulations [32t |33j. One of the most simple but nevertheless 
effective method of the bridge functional parameterization is to approximate the bridge functional 
by the exponential function. In the work [32] such kind of a parameterization was used to improve 
the predictions of the thermodynamical properties of organic compounds. In the recent work Q 
it was shown that by using the exponential bridge functional it is possible to predict the solute- 
solvent RDFs of the LJ balls of different sizes dissolved in the LJ fluid. The parameters of the 
fluid were chosen to coincide with the water parameters at normal conditions. This allows us to 
suggest that the exponential bridge functionals can be effective also for the aqueous solutions of 
LJ spheres. Moreover, in classical molecular dynamics (MD) simulations the carbon nanomaterials 
are described by the set of balls with L J potential which allows us to expect the good performance 
of the exponential bridge functional also in this case. In the current work the effectiveness of the 
exponential bridge functional for the aqueous solutions of carbon nanomaterials is tested. To fit 
the parameters of the exponential bridge function the infinitely diluted aqueous solutions of the L J 
spheres of different size are used. The bridges are parameterized in oder to have the best coincidence 
of the solvent DDF which are obtained from the 3DRISM equations with the solvent DDF which 
are obtained from the MD simulations. 



2. Method 

2.1 . Exponential Bridge functional for the 3DRISM equations 

In the current paper the 3DRISM model is used in oder to describe the infinite diluted solu- 
tions of LJ balls and carbon nanomaterials. The investigated systems contain one solute molecule 
surrounded by the solvent. The solute molecule is the three-dimensional object while the solvent 
molecules are represented in the RISM approximation as the set of interacting sites. 

The 3DRISM equations are written in the following way: 

iV 8 jte f 

Mr) = £ J c s (r-r')x« a (|r'|)dr' (2.1) 

£ =1 R3 

where h a (r), c a (r) are the total and the direct correlation functions of the solvent site a , X£a( r ) i s 
the susceptibility function of the solvent sites £ and a. The susceptibility function x?a( r ) i s defined 
in the following way: 

Xa(,(r) = u a(L (r) + ph ai (r) (2.2) 

where p is the solvent number density, h a ^(r) is the total correlation function of solvent sites £ and 
a, Lu a (: is the solvent intramolecular function which defines the structure of the solvent molecule. 
The solvent intramolecular function for the rigid molecules is defined in the following way: 

u (a (r) = 6 (a 5(r) + (I - 6 ia ) 6{r ~ r 2 * a) (2.3) 

where S^ a is the Kronecker's delta, S(r) is the Dirac delta function, r^ a is the distance between the 
sites £ and a in the solvent molecule. 

The 3DRISM equations (|2.1[) are completed with the closure relation: 

Mr) + 1 = e -^(r+Mr)- C «(r)+B«(r) ^ 

where U a (r) is the potential of the interaction of the solute with the solvent site a, B a (r) is the 
bridge functional of the solvent site a. 
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As it was mentioned above the bridge functional in general case is practically incomputable. 
In the paper [9] for the solutions of the LJ balls in the LJ fluid the following exponential bridge 
functional was used: 

B(r) = - exp(fc(r - r ) - CV so i ute /crsoivent) (2.5) 

where fc=30 nm~l, C=1.18, cr so i u t c and <7 so ivent & re the a LJ parameters of solute and solvent atoms 
correspondingly, ro is the point where U(ro) — 13.8fcsT. In the current work we use analogous 
bridge functional for the LJ spheres' aqueous solutions. Due to the spherical symmetry the functions 
in the equations (|2.1[) . (|2.4p depend only on the distance to the solute's center r = |r|. In the 
current work we use the following bridge functionals: 

B a {r) = - exp(fc(r - r£) - Ca/a Q ) (2.6) 

where oo = 0-316 nm is the oxygen a parameter in the SPC/E water model. 

2.2. Investigated systems 

In the current paper two classes of systems are investigated : (i) infinitely diluted aqueous 
solutions of the Lennard- Jones (LJ) spheres and (ii) infinitely diluted aqueous solution of the 
carbon nanotubes(CNT). The a LJ parameters of the spheres was a = Kao, where K={ 0.1, 0.2, 
0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.3, 1.4, 1.5, 1.8, 2 } , a = 0.316 nm. For the investigated 
systems the RISM/3DRISM calculations as well as the MD simulations were performed. For the 
MD simulations the program package Gromacs 4.5.3 was used[34]. LJ balls were put to the cubic 
box of size 4x4x4 nm filled with the water molecules. The SPC-E water model was used in the 
simulations [35j]. The simulations were performed in the NVT ensemble with the temperature 
T = 300K. 1 000 000 simulation steps of leap-frog integration method were performed. Each of the 
step corresponded to 1 femtosecond. For the RISM calculations the multi-grid algorithm described 
in the paper [5j was used. 

The second class of systems are the aqueous solution of carbon nanotubes (CNT). In that 
class investigate the single system was investigated: aqueous solution of the single-wall CNT. The 
nanotube of chirality (7,0) was used which corresponds to the diameter of 0.55 nm. The choice 
of the nanotube is explained by the fact that for the CNT of such chirality MD simulation data 
is available in the literature (see supporting information for works [36j, |37|). To test the effective- 
ness of the bridge functionals the 3DRISM calculations for the same system were performed. The 
LJ parameters of the OPLS force field was assigned to the carbon atoms of the CNT: gp=0.355 
nm, ec=0.07 kcal/mol. For the calculations the multi-grid 3DRISM algorithm was used [38]. The 
calculations with the exponential bridges as well as without bridges (for comparison) were per- 
formed. Because the 3DRISM equations without bridge cannot be solved numerically because of 
the divergence of the algorithm for the system without bridge the Kovalenko-Hirata closure was 
used [39]. This closure is the standard closure for the 3DRISM calculations (3^. In the 3DRISM 
calculations the Cartesian grid is used. It is characterized by two values: the buffer, which is the 
minimal distance from the solute atoms to the boundaries of the grid and the spacing which is the 
minimal distance between the grid points. For the 3DRISM calculations the grid with the buffer 
1.5 nm and spacing 0.02 nm was used. In the RISM and 3DRISM calculations the MSPC/E water 
model was used [40}. The water hydrogen had the following LJ parameters: o\fl-=0.8A , e#=0.046 
kcal/mol. 

2.3. Parameterization of the bridge functional 

In the work Q the optimal values for k and C parameters in the equation (|2.5p were determined. 
Because of the similarity in the solvent L J parameters in the current paper the same values as in the 
work [gj are used: fc=30 nm~l, C=1.18. The t-q parameter depends on the shape of the potential 
and LJ parameters of the atom pair, and thus is not universal. In the current work the optimal 
values of fg are determined and connected to the universal parameters, such as the positions 
of minimum of the potentials Uq{t), Ufj(r) which describe the interaction of the water oxygen 
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and water hydrogen with the solute. For the parameterization is more convenient to use not the 
functions Uq{t), Ujj(r) but dimensionless functions exp(— PUq{t)), exp(— /3U^(r)). They are the 
first approximations of the RDF functions (if one considers only the direct atom-atom interactions 
neglecting the influence of the solvent). The maximums' positions of the functions exp(— /?[/,§ (r)) 
and exp(— /3Ufj(r)) in the current work are denoted as r^Q and r a jj correspondingly. In the current 
work the correlation of the parameters r$, r$ with the r^Q, r c *£ for different sizes of the solute is 
investigated. 

It is supposed that the following relations hold: 

r$=ACS,r»=AC£ (2.7) 

where the constant A is the same for the water oxygen and water hydrogen, but can be different 
for different sizes of the solvent LJ ball: A = A(a so i uie /a so i ven t) 

For each solute size <7 so i u te the values of A from A — 0.5 to A = 1 are tested. For each of the A 
values the quality of the local solvent structure prediction is estimated. The center of mass of the 
water molecule nearly coincides with the oxygen atom which allows to estimate the local solvent 
structure prediction using the water oxygen RDF functions. The position and the height of the first 
peak of the water oxygen RDF as the numerical characteristic of the prediction of the RDF were 
used. These two parameters are very important for the thermodynamics of solution. It is reasonable 
to suggest that among these two characteristics the peak's position is more important. It describes 
the structure of the solvation shell while the peak's height is only the quantitative characteristic 
of the solvation shell's density. Thus, the first peak's position was used in order to find the best 
r o 'i r o values. After that it was tested how good the obtained bridge functional predicts the first 
peak's height of the RDF. 

3. Results 

RISM calculations for the infinitely diluted aqueous solutions of LJ sphere were performed. In 
the calculations the bridge functionals (|2.5|) were used. The rff, Tq were chosen using the formula 
(12.71) and the parameter A was varying in the range from A — 0.5 to A = 1. 

In the Figure Q] the dependence of the first peak position of the water oxygen RDF on the ratio 
^ = r o / r aO ^ s presented. The zones where the first peak position of the RDF obtained from the 
RISM coincides with the first peak position obtained from the MD simulations are marked with 
the gray color. As it can be seen, for the LJ spheres of small size the first peak of the water oxygen 
RDF weakly depends on A while for larger spheres the optimal value of A slightly fluctuates around 
the mean value. This allows to suggest that the approximation A = const can quite good fit to the 
MD-simulation's results. To determine the optimal value of A for each value of A the total error 
as a sum of the errors for different solute sizes was found: 

error{A) = ^ \rf D - r* ISM (A)\ (3.1) 

where r^ ID is the first peak position of water oxygen RDF obtained from the MD simulations, 
r RiSM j g first peak position of water oxygen RDF obtained from the RISM equations (|2.1j) 
and the closure relation (12. 41) . where the bridge functionals are calculated by using the formula 
(|2.6[) . and the parameters r® (A), r$ (A) are connected to the positions of the first maximum of 
the functions exp(— flU^r)), exp(— f3U^(r)) by the expressions (|2.7|) . 

In the Figure [5] the dependence of the total error of the first peak position prediction (|3.ip 
on the parameter A is presented. The minimal error corresponds to the value A = 0.77. So, the 
optimal universal bridge functional is the following: 

B aa (r) = - cxp(fc(r - A ■ r e *») - Ca/a ) (3.2) 

where A = 0.77, r°* p is the position of maximum of the function exp(— f3U^(r)). 

To test the bridge functional (|3.2[) the RISM calculations of the infinitely diluted aqueous 
solutions of LJ spheres of different size with the bridge functional (13.211 and without the bridge 
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Figure 1. Dependence of the first peak position of the water oxygen RDF calculated from the 
RISM equations with the exponential bridge (|2.6|l on the ratio A = r° jr^Q. The zones where 
the first peak position of the RDF obtained from the RISM coincides with the first peak position 
obtained from the MD simulations are marked with the gray color. 



1 




Figure 2. Dependence of the total error of the first peak's position of the water oxygen RDF 
prediction (pTTj) on the ratio A = r° /r e Jg. 
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Figure 3. The the water oxygen RDFs obtained from the RISM equations with bridge, from 
the RISM equations without bridge and from the MD-simulations. Three groups of RDFs for 
spherical LJ solutes of different sizes are presented: a so i uto =0.5cro, f so i uto =(Jo, <7 so i utc =1.5<7o, 
where ctq = 0.316 nm is the a LJ parameter of the water oxygen in the SPC/E water model. 
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1.820 
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0.436 
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0.270 


2.0 


1.794 
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0.414 
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0.306 
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1.876 


2.229 


0.352 


2.009 
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Table 1 . Comparison of the first peak's height of the oxygen water RDFs obtained by the three 
methods: the RISM simulations with bridge, the RISM simulations without bridge and the MD 
simulations. 
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0.00 
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0.49 
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0.49 


0.00 
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-0.021 




-0.003 



Table 2. Comparison of the first peak's position of the oxygen water RDFs obtained by the three 
methods: the RISM simulations with bridge, the RISM simulations without bridge, and the MD 
simulations. 



(B aa (r) = 0) were performed. The RDFs were compared to the RDFs obtained from the MD 
simulations. In Table [T] the comparison of the first peak's height is presented. In Table [5] the 
comparison of the first peak's position is presented. For the three selected sizes of the solute LJ 
spheres the oxygen water RDFs obtained from the RISM equations (with and without bridge) 
and from the MD simulations are presented in Figure [31 It can be seen that by using the bridge 
functional Q3.2p it is possible to predict the water oxygen RDF's first peak's position with the 
average accuracy of 0.003 nm (0.03A ). This is 3 times less than the resolution of the grid used 
for discretization the RDFs (O.lA ). Without using the bridge functional the average error of the 
peak's position prediction is 0.021 nm (0.2lA ) which is essential e.g. for the solvation free energy 
calculations. Introducing the bridge functional (I3.2[) also improves the predictions of the first peak's 
height of the oxygen water RDF. Without the bridge functional the peak's height is overestimated 
in average by 0.352 units (18.76%) while with the bridge functional the average error is only 0.133 
units (7.1%). 

In conclusion: introducing the bridge functional improves the accuracy of RISM calculations for 
aqueous solutions of the LJ spheres of different size. This also allows to suggest the effectiveness 
of the proposed bridge functional for the aqueous solutions of the carbon nanomaterials. To check 
this assumption the 3DRISM calculations for the aqueous CNT solution were performed with the 
proposed bridge functional and without it (using the Kovalenko-Hirata closure [39]). The bridge 
functional for the nanotube was calculated as the superposition of the bridge functionals of the 
atoms of the nanotube: 

B CK (r) = ^ J B CTQ (|r-r s |) (3.3) 

s 

where r s is the position of the atom s of the CNT, B aa (r) is defined by the formula (|3.2[) where 
a = (7(7=0.355 nm, and the sum is calculated over all the atoms of the nanotube. I calculated the 
mean number densities of the water oxygen and water hydrogen at the different distances from 
the nanotube axis. Comparison of the RISM calculations with the results of the MD-simulations 
which were done in the work [36j are presented in the Figure [U 

As it can be seen, introducing the bridge functional allows to predict almost accurately the 
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0.3 0.4 0.5 0.6 0.7 0.8 0.9 0.3 0.4 0.5 0.6 0.7 0.8 0.9 

distance from the CNT axis [nm] distance from the CNT axis [nm] 

(a) Water oxygen atoms' density distribution around (b) Water hydrogen atoms' density distribution 
the nanotube's axis around the nanotube's axis 

Figure 4. Comparison of the RISM calculations (with the exponential bridge and without the 
bridge) with the results of the MD-simulations for the CNT aqueous solution 

distribution of the water oxygen around the axis of CNT. The center of mass of the water molecule 
nearly coincide with the center of oxygen atom. That means that the 3DRISM with the proposed 
bridge functional can also correctly predict the water density around the carbon nanotube. The 
prediction of the distribution of the hydrogen atoms is not perfect. This can be explained by 
the fact that the 3DRISM describes the structure of the solvent molecule not enough accurate. 
However, it is necessary to notice that in the vicinity of the nanotube (at the distances <0.6 nm 
from the nanotube's axis) the water hydrogen distribution functions obtained from the 3DRISM 
calculations nearly coincide with the water hydrogen distribution functions obtained from the MD 
simulations. This fact allows to conclude that the bridge functional allows to predict correctly the 
water hydrogen interaction with the CNT. 

4. Conclusions 

In the current work the exponential bridge functional in the closure of the 3DRISM equation was 
parameterized and tested. Two classes of systems were investigated: the infinitely diluted aqueous 
solution of the LJ spheres and the infinitely diluted aqueous solution of the carbon nanotubes. 
The following a parameters of the LJ spheres were used: a = Kao, where K = { 0.1, 0.2, 0.3, 
0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.2, 1.3, 1.4, 1.5, 1.8, 2 }, cr = 0.316 nm. For all of the systems 
the MD simulations and the RISM calculations with the bridge functional (|2.6p and different 
values of Tq were performed. It was shown that introducing the bridge functional can essentially 
improve the prediction of the the water oxygen RDF's first peak's position. The dependency of the 
optimal values of Tq on the size of the LJ solute was investigated. It was shown that t-q is related 
to the position of the maximum of the function exp(— /3U aa (r)), where U aa (r) is the interaction 
potential of the LJ sphere of size a with the water's atom a. The optimal value of constant A 
in the relation (|2.7p was determined. As a result the universal functional (|3.2p was proposed. It 
was shown that by using the bridge functional Q3.2J) it is possible to predict the the water oxygen 
RDF's first peak's position with the accuracy of 0.003 nm (0.03A ). This is 3 times less than the 
resolution of the grid which was used for discretization of the RDFs. The RISM equations without 
the bridge give a systematic error of the first peak's position prediction of 0.021 nm (0.2lA ). 
It was shown that using the bridge functional Q3.2p also improves the predictions of the first 
peak's height of the water oxygen RDFs. Without the bridge functional the first peak's height is 
overestimated in average by 18.76% while with the bridge functional the average error is only 7.1%. 
The effectiveness of the proposed bridge functional was also tested on the more complex system: 
the aqueous carbon nanotube solution. For the investigations the carbon nanotube of the chirality 
(7,0) (diameter 0.55 nm) was used. The water oxygen and water hydrogen density distributions 
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at the different distances from the carbon nanotube were calculated. It was shown that the bridge 
functional allows almost accurately predict the water oxygen atom's distribution around the CNT 
(and thus the water molecules' distribution, because the center of mass of the water molecule is 
nearly coincides with the center of the oxygen atom). It was shown that the bridge functional 
allows to obtain the correct water hydrogen distribution in the vicinity of the CNT. 
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MicTKOBa cf)yHKi4m TpMBMMipHoV MOfle/ii B3aeMOflitOHMx 1461-iTpiB 
(3DRISM) ftnn BOflHux po3HMHiB Bynnei4eBi/ix HaHOMaTepia/iiB 

B.n. CeprieBCbKi/ifi 

Max Planck Institute for Mathematics in the Sciences, 22 Inselstrasse, 04103 Leipzig, Germany 

y CTaTTi 3anpOnaHOBaHO HOBwK MicTKOBl/lirt CpyHKi4i0Ha.il p/IX piBHflHb 33 M KHeHHfl TpHBUMipHOI MOfle/li B3ae- 
MOflitOHWX LieHTpiB (3DRISM). TeCTyCTbCiR eCpCpeKTH BHicTb MicTKOBOrO CpyHKL4i0Ha/iy BOflHWX p03HI/lHiB 

Byr/ieLieBwx HaHOMaTepiajiiB. y CTa~rri po3r.nflflaioTbCfl flBa K/iaca cwcTewi: (1) HecKiHHeH ho po36aB/ieHi BOflHi 
po3HMHM JleHHapfl-fl>KOHCOBCbKwx ccpep (2) HecKi h HeHHO po36aB.neHi BOflHi po3HMHM ByrneLieBMX HaHOTpy- 
6ok. MicTKOBMM cpy H KujoHa/i na pa MeTepi/i3ye~rbCfl Buxoflam/i i3 flaHHWx cwMyjiaujti. rioKa3aHO, mo fl/ia ycix 

flOCflifl>KeHWX CMCTeM M icTKOBM M cpyHKLHOHa/i MO>Ke 6yTH Ha6/lH>KeHWW eKCnOHeHL4iajlbHOK) cpyHKLiieto, mp 
3ajie>KWTb Ti/lbKM Bifl Cni BBiflHOLUeHHS p03MipiB p03 H l/l H e H Ol' peHOBWHW Ta p03H M HHM K3 . l~loKa3aHO, U4O BM- 

KopucTOBytoH n 3a nponoHOBa h ia v\ MicTKOBi/iK cpyHKL4ioHa/i e mo>k/imbhm: (1) tohho 3aB6an mth nojio>KeHHa 
nepiuoro niKy cpyHKujw po3nofli.neHHn mmibHOCTi (<t>PLU,) ki/ichio BOflw HaBKO.no po3HWHeHoT peHOBWHW (2) 

nOKpaiUMTW TOHHicTb 3aB6aHeHb BMCOTH nepLUOTO niKy 0PU4 KHCHKD BOflM HaBKO/lO p03HMHeH0T peHOBMHH 

(3) npaBM/ibHO 3aB6aHWTM noBefliHKy OPLH BOflHio BOfli/i no6iii/i3y Ha HOTpy6KM . 

K/iiOHOBi c/ioBa: 3DRISM, mIctkobmm dpyHKn'ionan , Byrneu,eBi HaHorpy6KH, dpyHKuJi po3nop,\neHHft 

lHIAJlbHOCTI 
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